################################################################################
# Libraries and Import
################################################################################
rm(list=ls())

# Libraries
library(sandwich)
library(tidyverse)
library(ri2)
library(stargazer)

# Import
load(file="data/dat_analysis_CN.RData")
cn <- dat_vignette

# Set seed
set.seed(30) # Randomization inference values will change based on seed

# Control number of simulations in randomization procedures
sims = 10000

################################################################################
# Define heterogenous groups
################################################################################
# Redefine education variable to binary college or no college
cn$edu <- with(cn, ifelse((edu == 1 | edu == 2), 0, 1))

# Liberal peace variables
cn$Japan_trade <- cn$Japan_export + cn$Japan_import
cn$inter_pos <- as.numeric(cn$international_sys==4) + 
  as.numeric(cn$international_sys==5)

# Create age and age-squared variables
cn$age <- as.numeric(as.character(cn$age_group))
cn$age2 <- cn$age^2

# Nationalism
cn$perception_fp = as.numeric(as.character(cn$perception_fp))
cn$international_sys = as.numeric(as.character(cn$international_sys))
cn$territory = as.numeric(as.character(cn$territory))
cn$China_citizen = as.numeric(as.character(cn$China_citizen))
cn$critique = as.numeric(as.character(cn$critique))

dv.names <- c("perception_fp", "territory", "China_citizen", "critique")
dv.names <- c("territory", "China_citizen", "critique")

factor.obj <- princomp(cn[, dv.names], cor=TRUE)
cn$nat_pca <- as.vector(factor.obj$scores[,1])

################################################################################
# Split each country sample into two datasets by treatment group
################################################################################
# Convert China response variable from factor to numeric
cn$D1 = as.numeric(as.character(cn$D1))

# Split dataset by treatment group - China
nat = cn %>% filter(vignette_group == 1 | vignette_group == 2)
econ = cn %>% filter(vignette_group == 3 | vignette_group == 4)

nat$z = with(nat, ifelse(vignette_group == 1, 0, 1))
econ$z = with(econ, ifelse(vignette_group == 3, 0, 1))

################################################################################
# Standard interaction tests
################################################################################
# Nationalist treatment
nat_party = lm(D1 ~ z + party + z:party, data = nat)
nat_education = lm(D1 ~ z + edu + z:edu, data = nat) 
nat_trade = lm(D1 ~ z + Japan_trade + z:Japan_trade, data = nat)
nat_nationalist = lm(D1 ~ z + nat_pca + z:nat_pca, data = nat)

# Economic treatment
econ_party = lm(D1 ~ z + party + z:party, data = econ) 
econ_education = lm(D1 ~ z + edu + z:edu, data = econ) 
econ_trade = lm(D1 ~ z + Japan_trade + z:Japan_trade, data = econ)

################################################################################
# Export tables with standard interaction test
################################################################################
# Nationalism 
stargazer(nat_nationalist, 
          out = "tables/TA7. china_nat_nationalist.tex", 
          se = starprep(nat_nationalist),
          title="Nationalism and support for nationalist interference (China)",
          label = "tab:nat_nationalist",
          align=TRUE, 
          dep.var.labels=c("Government Support"),
          keep = c("z", "nat_pca", "z:nat_pca", "Constant"),
          covariate.labels=c("Government interference", 
                             "Nationalist", 
                             "Government interference x nationalist", 
                             "Constant"),
          no.space=TRUE,
          keep.stat = c("n"))

# Liberal peace: education and trade
stargazer(nat_trade, nat_education, econ_trade, econ_education,
          out = "tables/TA12. china_edu_trade.tex",
          se = starprep(nat_trade, nat_education, econ_trade, econ_education),
          title="Testing the Liberal Peace Logic (China)",
          label = "tab:china_edu_trade",
          align=TRUE,
          dep.var.labels=c("Government Support"),
          column.labels=c("Nationalist", "Nationalist", "Economic", "Economic"),
          keep = c("z", "Japan_trade", "z:Japan_trade", "edu", "z:edu", "Constant"),
          covariate.labels=c("Government interference",
                             "Trade with Japan",
                             "Government interference x trade",
                             "College education",
                             "Government interference x college",
                             "Constant"),
          no.space=TRUE,
          keep.stat = c("n"))

################################################################################
# Individual CATEs method: Government action support by Party
################################################################################
# Nationalist: Calculate difference in CATEs by CCP support
nat_party_1 = lm(D1 ~ z, data = nat, subset = party == 1)
nat_party_2 = lm(D1 ~ z, data = nat, subset = party == 0)
dic <- abs(nat_party_1$coefficients[[2]] - nat_party_2$coefficients[[2]])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- nat
dic_null <- vector(mode = "integer", length = sims)

# 1. # Form full schedule of potential outcomes assuming constant effects
null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 
  
  # 2. Assign subjects randomly to treatment or control
  null$z_null = complete_ra(nrow(null))
  
  # 3. Calculate difference between estimated CATEs
  nat_party_null_1 = lm(D1_null ~ z_null, data = null, subset = party == 1)
  nat_party_null_2 = lm(D1_null ~ z_null, data = null, subset = party == 0)
  dic_null[[i]] <- abs(nat_party_null_1$coefficients[[2]] - nat_party_null_2$coefficients[[2]])
}

# 4. Calculate two sided p-value
p_nat_party <- sum(abs(dic_null) >= dic)/length(dic_null)

############# Economic: Calculate difference in CATEs by CCP support ###########
econ_party_1 = lm(D1 ~ z, data = econ, subset = party == 1)
econ_party_2 = lm(D1 ~ z, data = econ, subset = party == 0)
dic <- abs(econ_party_1$coefficients[2] - econ_party_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- econ
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 
  null$z_null = complete_ra(nrow(null))
  
  nat_party_null_1 = lm(D1_null ~ z_null, data = null, subset = party == 1)
  nat_party_null_2 = lm(D1_null ~ z_null, data = null, subset = party == 0)
  
  dic_null[[i]] <- abs(nat_party_null_1$coefficients[2] - nat_party_null_2$coefficients[2])
}

# Calculate two sided p-value
p_econ_party <- sum(abs(dic_null) >= dic)/length(dic_null)

################################################################################
# Government action support by education level
################################################################################
# Nationalist: Calculate difference in CATEs by education (college or no)
nat_edu_1 = lm(D1 ~ z, data = nat, subset = edu == 1)
nat_edu_2 = lm(D1 ~ z, data = nat, subset = edu == 0)
dic <- abs(nat_edu_1$coefficients[2] - nat_edu_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- nat
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 
  null$z_null = complete_ra(nrow(null))
  
  nat_edu_null_1 = lm(D1_null ~ z_null, data = null, subset = edu == 1)
  nat_edu_null_2 = lm(D1_null ~ z_null, data = null, subset = edu == 0)
  
  dic_null[[i]] <- abs(nat_edu_null_1$coefficients[2] - nat_edu_null_2$coefficients[2])
}

p_nat_edu <- sum(abs(dic_null) >= dic)/length(dic_null)

############# Economic: Calculate difference in CATEs by education #############
econ_edu_1 = lm(D1 ~ z, data = econ, subset = edu == 1)
econ_edu_2 = lm(D1 ~ z, data = econ, subset = edu == 0)
dic <- abs(econ_edu_1$coefficients[2] - econ_edu_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- econ
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 
  null$z_null = complete_ra(nrow(null))
  
  nat_edu_null_1 = lm(D1_null ~ z_null, data = null, subset = edu == 1)
  nat_edu_null_2 = lm(D1_null ~ z_null, data = null, subset = edu == 0)
  
  dic_null[[i]] <- abs(nat_edu_null_1$coefficients[2] - nat_edu_null_2$coefficients[2])
}

p_econ_edu <- sum(abs(dic_null) >= dic)/length(dic_null)

################################################################################
# Government action support by trade relationship
################################################################################
# Nationalist: Calculate difference in CATEs
nat_trade_1 = lm(D1 ~ z, data = nat, subset = Japan_trade == 1)
nat_trade_2 = lm(D1 ~ z, data = nat, subset = Japan_trade == 0)
dic <- abs(nat_trade_1$coefficients[2] - nat_trade_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- nat
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 
  null$z_null = complete_ra(nrow(null))
  
  nat_trade_null_1 = lm(D1_null ~ z_null, data = null, subset = Japan_trade == 1)
  nat_trade_null_2 = lm(D1_null ~ z_null, data = null, subset = Japan_trade == 0)
  
  dic_null[[i]] <- abs(nat_trade_null_1$coefficients[2] - nat_trade_null_2$coefficients[2])
}

p_nat_trade <- sum(abs(dic_null) >= dic)/length(dic_null)

############# Economic: Calculate difference in CATEs by education #############
econ_trade_1 = lm(D1 ~ z, data = econ, subset = Japan_trade == 1)
econ_trade_2 = lm(D1 ~ z, data = econ, subset = Japan_trade == 0)
dic <- abs(econ_trade_1$coefficients[2] - econ_trade_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- econ
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 
  null$z_null = complete_ra(nrow(null))
  
  nat_trade_null_1 = lm(D1_null ~ z_null, data = null, subset = Japan_trade == 1)
  nat_trade_null_2 = lm(D1_null ~ z_null, data = null, subset = Japan_trade == 0)
  
  dic_null[[i]] <- abs(nat_trade_null_1$coefficients[2] - nat_trade_null_2$coefficients[2])
}

p_econ_trade <- sum(abs(dic_null) >= dic)/length(dic_null)

################################################################################
# Government action support by pre-existing nationalism
################################################################################
# Nationalist: Calculate difference in CATEs by nationalism
nat_nat_1 = lm(D1 ~ z, data = nat, subset = nat_pca > mean(nat_pca))
nat_nat_2 = lm(D1 ~ z, data = nat, subset = nat_pca < mean(nat_pca))
dic <- abs(nat_nat_1$coefficients[2] - nat_nat_2$coefficients[2])

# Use randomization inference to test null that CATEs in both groups equal ATE
null <- nat
dic_null <- vector(mode = "integer", length = sims)

null$ate = lm(D1 ~ z, data = null)$coefficients[2]
null$D1_null = with(null, ifelse(z == 1, D1, D1 + ate))

set.seed(30)
for (i in seq_along(dic_null)) { 
  null$z_null = complete_ra(nrow(null))
  
  nat_nat_null_1 = lm(D1_null ~ z_null, data = null, subset = nat_pca > mean(nat_pca))
  nat_nat_null_2 = lm(D1_null ~ z_null, data = null, subset = nat_pca > mean(nat_pca))
  
  dic_null[[i]] <- abs(nat_nat_null_1$coefficients[2] - nat_nat_null_2$coefficients[2])
}

p_nat_nat <- sum(abs(dic_null) >= dic)/length(dic_null)

################################################################################
# Calculate adjusted p values
################################################################################
options(scipen=999)

Treatment = c("Nationalist", "Nationalist", "Nationalist", "Nationalist", 
              "Economic", "Economic", "Economic")

Covariate = c("Party", "Education", "Trade", "Nationalism",
              "Party", "Education", "Trade")

Unadjusted = c(p_nat_party, p_nat_edu, p_nat_trade, p_nat_nat, p_econ_party, 
               p_econ_edu, p_econ_trade)

Unadjusted = round(Unadjusted, 4)

Bonferroni = round(p.adjust(Unadjusted, method = "bonferroni"), 4)

BH = round(p.adjust(Unadjusted, method = "BH"), 4)

Holm = round(p.adjust(Unadjusted, method = "holm"), 4)

corrections = cbind(Treatment, Covariate, Unadjusted, BH, Holm, Bonferroni)

# Output table
stargazer(corrections, 
          header=FALSE, 
          summary = FALSE,
          column.sep.width = "2",
          title = 'Multiple comparisons corrections China',
          label = "corrections_china",
          out = "tables/TA8. corrections_china.tex")

################################################################################
# EXPLORATORY ONLY: Understanding Party Affiliates’ Tolerance
################################################################################
nat_critique = lm(D1 ~ z + party + z:party + critique + z:critique, data = nat)
nat_critique_cov = lm(D1 ~ z + party + z:party + prov + age + age2 + 
                              male + edu + critique + z:critique, data = nat)

nat$hardcore = ifelse(nat$China_citizen == 5, 1, 0)
nat_nationalist = lm(D1 ~ z + party + z:party + hardcore + z:hardcore, data = nat)
nat_nationalist_cov = lm(D1 ~ z + party + z:party + prov + age + age2 + 
                                 male + edu + hardcore + z:hardcore, data = nat)

stargazer(nat_critique, nat_critique_cov, nat_nationalist, nat_nationalist_cov,
          out = "tables/TA11. understanding_party_affiliates.tex", 
          title="Understanding Party Affiliates’ Tolerance of Government Interference with Protests (Chinese sample)",
          label = "understanding",
          align=TRUE, 
          dep.var.labels=c("Government Support", "Government Support", 
                           "Government Support", "Government Support"),
          keep = c("z", "party", "z:party", "critique", "z:critique", 
                   "hardcore", "z:hardcore", "Constant"),
          covariate.labels=c("Government interference", 
                             "CCP party affiliate", 
                             "Political conservatism",
                             "Hardcore nationalist",
                             "Government interference x Party affiliate", 
                             "Government interference x Political conservatism",
                             "Government interference x Hardcore nationalist",
                             "Constant"),
          no.space=TRUE,
          keep.stat = c("n","rsq", "adj.rsq"),
          add.lines = list(c("Covariate Adjustment", "No", "Yes", "No", "Yes")))

################################################################################
# EXPLORATORY ONLY: Liberal peace theory exploratory analysis
################################################################################
nat_model_trade = lm(D1 ~ z + Japan_trade + z:Japan_trade, data = nat)
nat_model_trade_cov = lm(D1 ~ z + prov + age + age2 + male + edu + Japan_trade + z:Japan_trade, data = nat)

nat_model_influence = lm(D1 ~ z + Japan_influence + z:Japan_influence, data = nat)
nat_model_influence_cov = lm(D1 ~ z + prov + age + age2 + male + edu + Japan_influence + z:Japan_influence, data = nat)

nat_model_inter = lm(D1 ~ z + inter_pos + z:inter_pos, data = nat)
nat_model_inter_cov = lm(D1 ~ z + prov + age + age2 + male + edu + inter_pos + z:inter_pos, data = nat)

# Output table using stargazer
stargazer(nat_model_trade_cov, nat_model_influence_cov, nat_model_inter_cov,
          out = "tables/TA14. china_liberal_peace.tex", 
          title="Testing the Liberal Peace Logic (China)",
          label = "libpeacechina",
          align=TRUE, 
          dep.var.labels=c("Government Support", "Government Support", "Government Support"),
          keep = c("z", "Japan_trade", "z:Japan_trade", "Japan_influence", 
                   "z:Japan_influence", "inter_pos", "z:inter_pos", "Constant"),
          covariate.labels=c("Government interference", 
                             "Trade with Japan", 
                             "Government interference x Trade with China", 
                             "Japanese products",
                             "Government interference x Japanese products",
                             "Internationalist",
                             "Government interference x Internationalist",
                             "Constant"),
          no.space=TRUE,
          keep.stat = c("n","rsq", "adj.rsq"),
          add.lines = list(c("Covariate Adjustment", "Yes", "Yes", "Yes")))